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, In this article we consider the linear stability of the spherically- 

■ symmetric stationary solutions of the Schrodinger-Newton equations. 

^ . These have been found numerically by Moroz et al |J and Bernstein 

et al |J. The ground state, characterised as the state of lowest en- 
ergy, turns out to be linearly stable, with only imaginary eigenvalues. 
The (n + l)-th state is linearly unstable having n quadruples of com- 
plex eigenvalues (as well as imaginary eigenvalues), where a quadruple 
consists of {A, A, —A, —A} for complex A. 

1 Introduction 

The Schrodinger-Newton equations (hereafter the SN equations) for a quantum- 
mechanical particle of mass m moving in its own gravitational potential, are 
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the pair of coupled non-linear partial differential equations: 



ih— = V^ + m$^ (1) 

at 2m 

V 2 $ = 4vrGm|^| 2 

where h is Planck's constant, \1/ is the wave function, $ is the gravitational 
potential, G is Newton's gravitational constant and t is time. 



The SN equations were introduced by Penrose ]l(J,[|ll[] in his discussion 
of quantum state reduction by gravity. Penrose suggested that a superposi- 
tion of two quantum states corresponding to two separated 'lumps' of matter 
should spontaneously reduce to one or other of the states within a time 
related to the self-energy of the gravitational field generated by the differ- 
ence of the densities. This idea requires that there be a preferred basis or 
special set of quantum states which collapse no further and these, Penrose 
suggested, should in the non-relativistic limit be the stationary states of the 
SN equations. 

The spherically-symmetric stationary states of the SN equations have 
been found numerically by a number of authors starting with Ruffini and 



Bonazzola |T^| and including Moroz, Penrose and Tod |J and Bernstein, 



Giladi and Jones 0. Our aim in this article is to test the linear stability 
of these solutions by linearising the time-dependent SN equations about a 
spherically-symmetric stationary state. In a later article, we shall present a 
numerical study of the full non-linear evolution in certain symmetric cases. 

The SN equations as introduced by Penrose are closely related to the 
Schrodinger-Poisson equations which have been studied for much longer (see 
e.g. 0). They are also the non-relativistic limit of the Einstein equations 
with a complex Klein-Gordon field as source (see e.g. fl2"f ) which is why their 
solutions are sometimes known as boson stars (see e.g. 



The plan of this article will be as follows. In the remainder of this section 
we review the nondimensionalised SN equations and collect some analytical 
results on them. In Section 2, we review the spherically-symmetric stationary 
states and in Section 3 we derive the linearised time- dependent SN equations. 
In Section 4 we describe the numerical methods which will be used to solve 
the linearised equations and in Section 5 we present the results. 

We shall find that the ground-state solution is linearly stable, having 
purely imaginary eigenvalues, and the n-th excited state, equivalently the 
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(n + l)-th state, has n unstable modes corresponding to n quadruples of 
complex eigenvalues, together with purely imaginary eigenvalues. These re- 
sults are in line with the expectation one might have from similar nonlinear 



problems like those described in [16| and 



Following Moroz et al || we begin by introducing dimensionless variables 
for the system (JJ). We write the rescaled $ as ip, <fi, and we rescale space 
and time but write them using the same variables as before to obtain the 
nondimensional SN equations: 

^ = -VV + #, (2) 



The most important rescaling, which we note for later use, is that for the 
time variable, for which t new = ytoid with 

32vr 2 G 2 m 5 

7 = (3) 

Stationary solutions as usual take the form ip(x, y, z)e~ lEt and satisfy the 
nondimensional time-independent SN equations: 

Etj) = -V 2 V> + #, (4) 
VV = H 2 - (5) 

At this point it is worth collecting some analytical results on the SN equa- 
tions. 

• Existence and uniqueness for solutions of (0) is established by the fol- 
lowing simplified version of a theorem of Diner et al . 

Given x( x ) £ i/ 2 (R 3 ) with Z/2-norm equal to 1, the system 
@ has a unique, strong solution if)(x, t), global in time, with 
V>O,0) = x(x) and J |^| 2 = 1. 

Illner et al @ give precise regularity properties for the resulting solu- 
tion. 
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The equation can be obtained as a variational problem from the 
Lagrangian: 

8 = /(v^-V^ + ^M 2 ) (6) 

where it is understood that is the solution of (^), and the overbar 
denotes complex conjugation. (All integrals will be over R 3 unless 
indicated to the contrary.) Here E is a Lagrange multiplier for the 
normalisation constraint 

" Ai 2 = 1- (7) 



The quantity E of (|6]) is bounded below ([]H[). It is reasonable to 
suppose, though a proof has not appeared in the literature, that the 
infimum of 8 is attained and that the ground state is the lowest en- 
ergy spherically-symmetric solution found numerically by the authors 
mentioned above and proved to exist in ||. 

The system (|2]) preserves the normalisation ([?[) and the quantity £ of 
(P) even with time-dependent ip and <j). We shall therefore call S the 
conserved energy. Note that it is different from the energy eigenvalue 
E for a stationary state appearing in (f|). 

If we define the kinetic energy T and potential energy V in the obvious 
way by 

T= /|V^| 2 , V = [m 2 , (8) 



then £ = T + \V . It was shown in |T4| that, in a stationary state, 



14 1 
T = — E, V = -E, £ = -E. (9) 

3 ' 3 ' 3 V ' 

In particular, as one expects, the energy eigen-values are all negative. 



• There is an interesting observation on the dispersion of the wave- 
function due to Arriola and Solar JT|. Define the second moment Q 
by 

Q = /|x| 2 H 2 (io) 
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assuming that the initial data is such that this integral exists, then it 
follows from (^j that 

Q = J (8|Wf + 20|</f) = 88 - IV. (11) 

Now recall that, as a consequence of the maximum principle (see e.g. 
H) cf) is everywhere negative and therefore so is V. Thus if the con- 
served energy £ is positive, then the dispersion grows at least quadrat- 
ically in time. In this sense, a positive energy solution necessarily 
scatters. 



2 The spherically-symmetric stationary states 

For a spherically-symmetric stationary state, there is no loss of generality in 
assuming that the wave-function is real. The time-independent SN equations 
can therefore be written in the form 

(rV)rr = -ripU, (12) 
(rU) rr = -rip 2 , 

where we have introduced the variable U = E — <ft and used a subscript r for 
d/dr. 

The boundary conditions are that ip — > as r — > oo and ip r = = <p r at 
r = 0. If (i/),U,r) is a solution then so is (A 2 ?/>, X 2 U, A _1 r) for any nonzero 
A. This scaling can be used to impose the normalisation (|7|) retrospectively. 
To solve ([12]), Moroz et al H] use a shooting method. They set U(0) = 



1, which can be done with a suitable A leaving the solution not correctly 
normalised, then they choose values of ip(0) and integrate away from r = 0. 
The integration continues until the solution for ip diverges at some value of 
r. This value of r can be pushed out to larger distances and the solution 
may diverge to large positive or large negative values. By adjusting -^(0) it 
is possible to bracket a value at which i/j remains finite and then E can be 
identified from the limit of U at large distances. 

With better triggering to recognise the diverging solutions, this technique 
can be refined and it becomes feasible to obtain the first 50 energy levels 
in a reasonable time. The first 20 eigenvalues calculated with the above 
routine are given in Table |l|. The first 16 agree to this order with those 
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Number of zeros 


Energy Eigenvalue 





-0.163 


1 


-0.0308 


2 


-0.0125 


3 


-0.00675 


4 


-0.00421 


5 


-0.00287 


6 


-0.00209 


7 


-0.00158 


8 


-0.00124 


9 


-0.00100 


10 


-0.000823 


11 


-0.000689 


12 


-0.000585 


13 


-0.000503 


14 


-0.000437 


15 


-0.000384 


16 


-0.000339 


17 


-0.000302 


18 


-0.000271 


19 


-0.000244 


20 


-0.000221 



Table 1: The first 20 eigenvalues 

given by Bernstein et al 0. A log-log plot of the n-th energy E n against 
n is shown in figure [I| The slope is asymptotically very close to -2, which 
would be the exact value for the Hydrogen atom (as was previously noted by 
Bernstein et al ||). Once one has the eigenvalues, one can alternatively use 
a spectral method to find the eigenf unctions and corresponding potentials at 
the Chebyshev points. These will be needed in Sections 3 and 4 (where the 
definition of Chebyshev points is also given). 



6 




Figure 1: Log- log plot of the energy values against n 

3 The perturbation equations 

In this section we shall set up the linear stability problem, ready for numerical 
solution in Section 4. We look for a solution to (Q) which can be expanded 
in terms of a small parameter e as 

V>(x,t) = Vo + eV'i + ---, (13) 

0(x, t) = (f) Q + 601 + . . . , 

Substituting ( p~3D into the SN equations (0) and equating powers of e we 
obtain 

#ot + V 2 V>o = ^o^o (14) 
V 2 o = l^ol 2 

iipu + V 2 V>i = ^o0i + ^i0o (15) 
V 2 0i = VoV'i + ^o^i- 
We restrict to spherical symmetry and take 

^ = i? (r)e^, (16) 
0o = £-t/o(r), 
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where Rq and Uq, both real, are known from the work described in Section 
2. We seek solutions of the form 



Vi = l -R{r,t)e- iE \ (17) 
r 

where R is complex and <j) is real. With this choice of variables, (|T5|) becomes 

iR t + R rr = Ro<p-U R, (18) 
(f) rr = RqR + RqR. 

We look for a solution of the form 

R = (A + B)e xt + (A-B)e xt , (19) 

where we assume for now that A is not real and that A, B and W are functions 
only of r. Now we substitute into fll8|) and equate coefficients of e A * and e xt 
still assuming that A is not real. This gives 

W rr = 2R A, (20) 

A„ + U A - RqW = -iXB, 
B rr + U B = -iXA. 

If (A, B, W) is a solution of (p0|) with eigenvalue A then (A, — £>, W 7 ) is a solu- 
tion with eigenvalue —A. Thus eigenvalues come in quadruples (A, —A, A, —A) 
for complex A. If there exists a solution of (p0|) with real A then the real and 
imaginary parts of the system decouple. We may assume without loss of 
generality that A and W are real and B is pure imaginary and then there 



is once again a solution of the form of (|19|) to dig) . Thus ( P0|) gives the 
perturbation equations with no restriction on A. From the definitions ( |17| ) 
we need A, B and W to vanish at r = and to vanish asymptotically as r 
tends to infinity (or at the outer edge of the grid in a numerical calculation). 



We note that A = will be an eigenvalue of ( P0|) with A = W = 0, 
B = R . This is a trivial 'zero-mode' corresponding to an infinitesimal con- 
stant rotation in the phase of the unperturbed state. It will be found by the 
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numerical calculation in Section 5 and we shall discard it by hand. 



We can obtain a condition on A 2 as follows: multiply the second of ( |20| ) 
by A and the third by B and integrate to find with the aid of the first that 

roo _ roc j 

-iX ABdr = / (Uq\A\ 2 - \A r \ 2 - -\W r \ 2 )dr, (21) 
Jo Jo 2 

-iX ABdr = / (U \B\ 2 - \B r \ 2 )dr. 
Jo Jo 



The right-hand-sides in (^Tj) are both real, so that the quotient 

-iX / °° ABdr 



(22) 



iX^ ABdr ' 

is real provided it is well-defined. Hence we have a dichotomy: 

either X 2 is real or I ABdr = 0. (23) 
Jo 

This will provide a useful check on the numerics in Section 5. 

By a more sophisticated argument, which we shall relegate to an ap- 
pendix, we may obtain an upper bound on the real part of any eigenvalue 
in terms of the energy of the state being perturbed. For perturbations of a 
state with energy-eigenvalue E this is 



|K(A)|<— V^E. (24) 

This is an interesting analytic result which limits the rate at which an un- 
stable state can decay. It also will provide a useful check in Section 5. We 
now turn to the problem of solving the perturbation equations. 



4 Solving the perturbation equations 

We solve the perturbation equations by a spectral method (see e.g. [|i~5f ) 
which reduces the calculation to that of solving a matrix eigenvalue problem. 

We choose a range (0, L) for r, regarding r = L as being at infinity for the 
purpose of setting boundary conditions, and set x = 2-^ — 1. Next we choose 
a number N and approximate any function f(x) by the unique polynomial 
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p(x) of degree N or less such that p(x)~) = Vj. for k = 0, 1 . . . N , where 
v k — f{xk) and Xk = cos(kir/N) are the Chebyshev points. 

To approximate the derivative we introduce Wk by Wk = p\xk) for k = 
0,1 ... N. Then the differentiation matrix for polynomials of degree N, de- 
noted Dn, is defined by the equation 



\w N ) 



D 



N 



( v \ 
V v N J 



(25) 
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For an explicit formula for D^, see e.g. 

The second- derivative matrix is just D 2 N . The requirement that A, B and 
W be zero at the boundary is imposed by deleting the first and last rows and 
columns of the relevant differentiation matrix, in this case Dj^. This yields 
an (N — 1) x (N — 1) matrix conveniently denoted by D 

For the perturbation equations 
value equation 



N 



we therefore obtain the matrix eigen- 



where 



-2R 








n 2 

N 



Un 



-D 2 




iX 











( A 


-/ 







B 





/ 


: 


\ W 



(26) 



and 





( Mxi) \ 




A(x N -i) 




B(x 1 ) 






\W ) 


B(x N ^) 




W(xx) 




{ W(x N ^) ) 



Ro = diag(Ro(xi),Ro{x2),...,Ro{x N - 1 )) 
U = diag(Uo(x 1 ),Uo{x 2 ),...,U (x N ^i)). 



(27) 



10 



where diag indicates a diagonal matrix. 

One might worry that since the matrix on the right in ( p6j) is singular, 
accuracy in the eigenvalues was lost. By solving the equation for W in terms 
of A we can rewrite (^6]) as 

D% + U \ (A \ 

D% + U - 2R (D%)- 1 R J [ B ) 

However the difference between the eigenvalues computed by using (p8|) and 
by using (|26| ) turns out to be small and we can just use (|26|). 

5 Numerical results 

We consider first the results obtained by solving (p6|) about the ground state 
of the system @,(^). Recall for comparison that the eigenvalue for the 
ground state in non-dimensionalized units is —0.163. The conversion factor 
to dimensional units is given by ([3]). In figure ^| we plot those eigenval- 
ues obtained by solving ( p6[ ) which are closest to the origin, excluding the 
near-zero eigenvalue which does not correspond to a non-trivial perturba- 
tion. To compute these results, we used N = 60 and L = 150. (That the 
near-zero eigenvalue found does indeed correspond to the trivial zero-mode 
can be seen by plotting the corresponding eigenfunction. One finds that A, 
W are very small (O(10~ 7 )) as compared to B, and that B is close to R$, so 
that this solution indeed corresponds to the zero-mode foreseen in Section 3.) 

In figure || we plot all the eigenvalues obtained by solving (|26|) (note the 
difference in vertical scale from figure 0) so that the quadruple of complex 
eigenvalues appears as a pair) to show that up to these limits there are no 
eigenvalues other than imaginary ones. The first few eigenvalues obtained, 
now including the near-zero one, are presented in table |2|. The non-trivial 
eigenvalues are all imaginary and we may conclude that the ground state is 
linearly stable. 

To test convergence of the method we plot graphs of a calculated eigen- 
value against increasing iV and L. As an example we show the graph of the 
fifth eigenvalue 0.0765i as the value of iV increases in figure f| and as the 
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-0.08 — 
— 1 



* 



-0.8 -0.6 -0.4 



0.2 0.4 0.6 



Figure 2: The smallest eigenvalues of the perturbation about the ground 
state 



-1 -0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8 



Figure 3: All the computed eigenvalues of the perturbation about the ground 
state 
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±0.000000116 
±0.0341i 
±0.0603z 
±0.0688z 
±0.0731z 
±0.0765i 
±0.0810i 
±0.0867i 



Table 2: Eigenvalues of the perturbation about the ground state 



value of L increases in figure |JJ There is satisfactory convergence. 

Using the same method we compute the numerical solutions of the pertur- 
bation about the second state. This time we obtain some eigenvalues with 
nonzero real parts. In figure ^| we plot the eigenvalues nearest to the origin 
for the case where N = 60 and L = 150. In figure [7] we plot all the eigenval- 
ues obtained with a different vertical scale (so that the quadruple of complex 
eigenvalues appears as a pair) to see that up to these limits there are no 
other complex ones. 

From Section 3 we have two analytic conditions on the eigenvalues. The 
first ( p3[ ) relates to complex eigenvalues and we do now have some. By the 
result of (|23|) we expect that for these Jq°° ABdr = 0. To see whether this is 
the case (up to numerical error) we compute 

n = \f° ABdr \ (29) 

U}\A\'dr)tU}\B\*dr)l 

which we want to be much less than one. For L = 145 and N = 60 we 
display in table the calculated values of Q with the eigenvalues nearest the 
origin. As expected, for the eigenvalues with nonzero real part, Q is zero 
within numerical error and the result of ( p3[) is confirmed. 

For the third state, with N = 60 and L = 450, we present the eigenvalues 
in figure ||. We have also calculated the eigenvalues for the fourth state and 
a consistent pattern emerges: there are n-quadruples of complex eigenvalues 
for perturbations about the (n + l)-th bound state (equivalently for the n-th 
excited state); the rest of the eigenvalues come in purely imaginary pairs. 

With the computed eigenvalues for the higher states available we can 



check their consistency with the other result, (24), from Section 3. In ta 
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x 1 O 7 
1 8 | 



60 
N 



Figure 4: The change in the sample eigenvalue with increasing values of iV 
(L = 150) 




100 110 120 130 140 150 160 170 180 190 

L 



Figure 5: The change in the sample eigenvalue with increasing values of L 
(N = 60) 
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-0.04 1 

-1 .5 



Figure 6: The lowest eigenvalues of the perturbation about the second bound 
state 




Figure 7: All the computed eigenvalues of the perturbation about the second 
bound state 
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A 


Q 


0.003* 
0.00860i 
-0.00139 -O.OlOi 
0.00139 + 0.010i 
0.0153z 


0.235 

0.368 
3.53exp(-14) 
3.92exp(-14) 

0.894 



Table 3: Q for different eigenvalues of the perturbation about the second 
state 



ble ^ we compare the maximum of dt(X) for the (n + l)-th state against 
the bound fl24l) . The computed eigenvalues comply comfortably with this 
analytic bound. When we have the eigenvalues for a particular state, we 



State 


Computed maximum 3?(A) 


Bound from (fPJ) 


1 





0.00362 


2 


0.00139 


0.00158 


3 


0.000520 


0.00101 


4 


0.000225 


0.000738 


5 


0.000114 


0.000583 


6 


0.0000653 


0.000482 



Table 4: Bound on the real part of the eigenvalues 

can check the accuracy of the spectral method by solving (|20| ) with each 
eigenvalue in turn. Using a Runge-Kutta method, we arrive at the same 
eigenfunctions as by the spectral method. 



To summarise, in this paper we have used a spectral method to analyse the 
linear stability of the stationary solutions of the SN equations in the context 
of the time-dependent SN equations. We have found that the ground-state 
solution is linearly stable having purely imaginary eigenvalues and the n-th 
excited state has n unstable modes corresponding to n quadruples of complex 
eigenvalues, together with purely imaginary eigenvalues. These results are 
tested for convergence and compliance with two analytical results. They 
are also in line with the expectation one might have from similar nonlinear 
problems like those described in [Tj| and Q. We shall consider the nonlinear 
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Figure 8: The eigenvalues of the perturbation about the third bound state 

stability in a second paper where this picture is confirmed and one can see 
the excited states decaying to the ground state. 
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Appendix 



In this appendix, we shall prove (plj ) starting from (|20|) . The calculation 
uses the best constant for the Sobolev inequality in three-dimensions given 
by Aubin (0). To work in three-dimensions, we first define a = A/r,b = B jr 
and w = W/r. In these variables we can write (^) as 



V 2 w 

V 2 a + Uqcl — Rqw 
V 2 6 + U b 



2RqCl, 
—iXb, 
iXa. 



(30) 



We multiply the second of (|30| ) by a, the complex conjugate of the third by 
b, add and integrate by parts to find 



i / [Xaa + Xbb] 



Robw, 



(31) 



where, as usual, the integrals are over R 3 . Now we make repeated use of the 
Holder inequality and the Sobolev inequality in 3-dimensions: 



(/ F 6 )i <K{( |VF| 2 )5 



where K is the relevant Sobolev constant, namely 



K 



2f 

32 7T3 



From the first of ( p0|) and the Holder inequality we find 



Vwl 2 < 2 / \Roaw\ < 2(/ \R 



and using the Sobolev inequality in this gives 

(/ \Vw\ 2 ) l i <2K(J \R \ 3 )HJ 



\a\ 2 )^ 



a| 2 )i 



kl 6 ) 1 



(32) 



(33) 
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Next we calculate 

/ \RoH < (/ \R \ 3 )hf \b\ 2 )Hj M 6 )* 

< K{J \R \ 3 )HJ \b\ 2 )Hf l^wl 2 ) 1 * 

< 2K 2 (J\R f)l(J\a\ 2 )Hj\b\^ 

where the first step uses the Holder inequality, the second uses the Sobolev 
inequality and the third uses equation (|33|) . We shall put this together with 
(|3T|). First we normalise the perturbation so that 

J \a\ 2 = cos 2 9, J \b\ 2 = sin 2 9, (34) 

for some 9 and set A = jj, + iv. Now from and fl3"4|) we have 

|/i + wcos20| < K 2 sin29( J \R \ 3 )^ ■ 

from which it follows that 

\^\<K 2 {J\R^)l. 

We need an upper limit for the right-hand-side in this, which we obtain as 
follows: 

J\R f < (JrdHJr^ 

< |Vi? | 2 )* 

3 3 

= K*T* 

using the Holder and Sobolev inequalities again and the definition (|]) of T. 
Using (Q) for T in a stationary state and the definition fl3"2"| ) of K we finally 
arrive at the desired limit 

N ^ q^^- ( 35 ) 
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